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ABSTRACT 

We present high resolution simulations on the impact of ionizing radiation of massive O-stars on 
the surrounding turbulent interstellar medium (ISM). The simulations are performed with the newly 
developed software iVINE which combines ionization with smoothed particle hydrodynamics (SPH) 
and gravitational forces. We show that radiation from hot stars penetrates the ISM, efficiently heats 
cold low density gas and amplifies over-densities seeded by the initial turbulence. The formation of 
observed pillar-like structures in star forming regions (e.g. in M16) can be explained by this scenario. 
At the tip of the pillars gravitational collapse can be induced, eventually leading to the formation of 
low mass stars. Detailed analysis of the evolution of the turbulence spectra shows that UV-radiation 
of O-stars indeed provides an excellent mechanism to sustain and even drive turbulence in the parental 
molecular cloud. 

Subject headings: stars: formation — ISM: structure — turbulence — ultraviolet: ISM — methods: 
numerical 



1. INTRODUCTION 

Some of the most spectacular structures in the molec- 
ular ISM are observed in the vicinity of hot O/B-stars 
or associations, e.g. the Horsehead nebula (B33), the 
three pillars of creation in M16 and the Elephant trunk 
(BRC 37) in IC1396. For the pillars in M16 lSugitani et all 
(2002) find a head to tail structure with the denser head 
pointing towards the OB stars of NGC661. In addition, 
young stellar objects (YSOs) are present at the tips of the 
pillars . In the Horsehead nebula IWard-Thompson et al.l 
( 2006) report two core- like structures that might undergo 
subsequ ent gravitational c ollapse. Very recent observa- 
tions bv llkeda et al.l (|2008h report several YSOs close to 
the tip of BRC37. As a common feature these pillar- 
shaped nebulae point towards a source of ionizing radi- 
ation and show signs of present or future star formation 
at their tips. 

Up to now the precise physical processes leading to 
the formation of these structures are not fully under- 
stood. The morphologies suggest that feedback effects 
of UV-radiation and winds of massive stars play an im- 
portant role in the formation of the pillars. In addition, 
the radiation might have a strong impact on the overall 
evolution of the parental cloud. Furthermore, molecular 
clouds are observed to be highly turbulent structures. 
There is evidence that this turbulence can support the 
clouds against gravitational collapse and thereby control 
star formation. As hydrodynamic and MHD turbulence 
decays rather quickly, the only way to explain this high 
level of turbulence would be to drive the turbulence - ei- 
ther on large scales by i.e. supernova explosions or on 
small scales from within the cloud by stellar outflows, 
win ds or ionization (see e.g. lElmegreen k Scald [2004 
and lMac Low k Klessenll2004L for reviews). The possi- 
bility of ionization driven turbu lence has been indicate d 
by e.g. semi-analytic models of iKrumholz et all (|2006h . 
In this Letter we test the hypothesis using high reso- 
lution numerical simulations with the newly developed 



code iVINE (jGritschneder et al.ll2009l hereafter G08). 

On the theoretical si de progress has been made since 
Elmc green et "all (| 19951 ) first presented two-dimensional, 
grid-based simulations showing that the expansion of an 
HII region into the surrounding ISM can trigger star for- 
mation by sweeping up the cold material. This is called 
'collect and collapse'. Another proposed scenario is the 
'radiation driven implosion', where pree xisting density 
structures are driven into collap se (see e.g. lBertoldilll989l 
iKessel-Devnet k Burkertll2003l and G08). 

For the numerical treatment of radiation in simula- 
tions several codes have been developed (see llliev et aTl 
120061 and references therein). Recent applications for 
the treatme nt of ionizing radiatio n in grid based codes 
includ e e.g. iMellema et alj (|2006h and IKrumholz et ail 
(|2007f h In SPH-cod es imp lementations have been pre- 
sented bvlDale et all pOOl . lPawlik k Schavd (l2008h _and 
lAltav et all (j2008h . Simulations by iDaleeLaJj ((20071) 
show that ionizing radiation can slightly enhance the for- 
mation of cores in a globally unbound molecular cloud 
of 10 4 M Q . With their choice of initial conditions the 
positive feedback, the additional or faster formation of 
cores, outweighs the negative feedback, the disruption 
of cores. All these applications calculate the effect of 
a point source on the surrounding medium, thereby fo- 
cussing much more on the global effect of the ionization. 
However, neither the detailed morphology of the gas nor 
the impact of the ionizing radiation on the turbulence 
has been investigated so far. 

2. INITIAL CONDITIONS 

We set up a box of gas with sides 4pc long at a 
temperature of T co u = 10K and a mean number den- 
sity of n = 300cm~ 3 , which resembles a slightly denser 
part of a molecular cloud. The gas mass in the box 
is 474M which corresponds to 25 Jeans masses. To 
mimic initial turbulence we employ a supersonic tur- 
bulent velocity field (Mach 10) with a steep power-law 
E(k) oc fc~ 2 , where only the largest modes k — 1..4 are 
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populated initially. This setup is allowed to freely decay 
under the influence of isothermal hyd rodynamics sim- 
ulated with the tree/SP H-code VINE (|Wetzstein et all 
120081 iNelson et al.|[2008h . The individual particle time- 
steps in VINE are determined by using an accuracy pa- 
rameter of T acc = 1.0 and a Courant-Friedrichs-Lewy 
(CFL) tolerance parameter of tcfl = 0.3. We also use 
an additional time-step criterion based on the maximum 
allowed change of the smoothing length with an accuracy 
parameter of Th = 0.15. 

After fss lMyr a Kolmogorov-like power-law with 
E(k) oc k~s is well established on all resolvable scales. 
The velocities now correspond to Mach 5. This initial 
setup with the turbulent velocities is shown in Fig. [T] 
(top panel) , the corresponding power-spectrum in Fig. [2] 
(top panel). 

With these turbulent initial conditions we perform 
two simulations, one with and one without the inclu- 
sion of ionization. To account for the UV-radiation of 
a young massive star we use iVINE (G08), a new par- 
allel implementation of ionizing radiation in the tree- 
SPH code VINE. Here we assume plane-parallel infall of 
UV-radiation onto one transparent side of the simulated 
area, which enables us to perform simulations at yet un- 
matched high resolution. From the surface the radiation 
is propagated by a ray-shooting algorithm. The size of 
the rays is determined by the smoothing-length of the 
SPH-particles, i.e. the width a particle occupies. Along 
these rays the radiation is calculated. This provides us 
with an ionization degree r\ for each SPH-particle, which 
is then used to assign a new temperature to each particle 
by linear interpolation. 

T = Thot • V + T co m ■ (1 - J]), (1) 

where T co id = 10K is the initial temperature of the cold, 
unionized gas and Thot = 10 4 K is the aver age tempera- 
ture of the ionized gas (see e.g. IShml Il99lh . The gas is 
assumed to be atomic hydrogen. Both gas components 
are close to thermal equilibrium since the heating and 
cooling timescales are much shorter than the dynamical 
timescales. We treat the gas with an isothermal equa- 
tion of state (7 = 1) as for the density range in our sim- 
ulations heating and cooling should balance each othe r 
to approximate isothermality (see e.g. lScalo et al.lfT998h . 
However, in reality th e situation is more c ompli cated. 
Recent simulations by iGlover fc Mac Low! (|2007f ) indi- 
cate an equation of state of the thermal equilibrium gas 
which is softer than isothermal (7 = 0.7 — 0.8). For a de- 
tailed prescription of the iVINE-code along with several 
analytical test cases see G08. In the simulations pre- 
sented here the radiation was calculated on more than 
(60) 2 rays, with the additional inclusion of five levels of 
refinement, leading to a spatial resolution of 2 x 10 _3 pc 
in the radiation. The photon flux per unit time and area 
is set to Fl v — 5 x 10 9 photons cm~ 2 s _1 , allowing the 
radiation to penetrate the first 10% of the cloud immedi- 
ately. This corresponds to setting up our simu lation to be 
right at the border of the Stromgren-sphere (Stromgren 
H9m which can be immediately ionized by an O-star or 
association. 

The radiation is impinging from the negative x- 
direction. Hydrodynamics is calculated with periodic 
boundaries in the y- and z-direction. The boundary is 
assumed to be reflecting in the negative x-direction to 
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Fig. 1. — Evolution of the turbulent ISM under the influence 
of UV-radiation impinging form the left hand side. Color coded 
is the surface density projected along the z-direction. The time of 
the snapshot is increasing from top to bottom. 

represent conservation of flux towards the star, whereas 
in the positive x-direction the gas is allowed to stream 
away freely. Gravitational forces are calculated without 
boundaries. This is valid as the free-fall time of the whole 
simulated area is « 3Myr, which is much larger than 
the simulation time. To ensure a correct integration of 
all quantities we use the individual time-stepping-scheme 
of VINE with the same parameters as for the freely de- 
caying turbulence (see above). For the tree-based calcu- 
lation of gravitational forces we use a mult i-pole accep- 
tance criterion (MAC, ISpringel et al.ll2~001h with a tree 
accuracy parameter of 9 = 5 x 10~ 4 . The correct treat- 
ment of the ionization and the resulting acceleration of 
the particles is obtained by a modified CFL-condition as 
discussed in G08. The simulations are performed with 
2 x 10 6 gas particles on a SGI Altix 3700 Bx2 supercom- 
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puter. The entire calculation took approximately 100 
wall clock hours on 16 CPUs. 

3. RESULTS 

3.1. Morphology and Formation of Cores 

At the beginning of the simulation the R-type front 
immediately reaches into the first 10% of the box, with 
the radiation penetrating further into the low density 
parts of the cold gas. After a hydrodynamical crossing 
timescale of the hot gas (thot ~ SOkyr) 1 the ionized gas 
reacts to its increase in temperature and starts to exert 
pressure on the cold gas. The cold gas is compressed and 
pushed away from the source, leading to a systematic ve- 
locity in the x-direction. At the same time the radiation 
has penetrated and ionized the ISM along channels of 
low density gas. Now these low density regions expand 
and start compressing the denser, unionized regions es- 
pecially tangential to the direction of radiation. Thus 
the preexisting density structures, which are seeded by 
the turbulent initial conditions get enhanced as shown in 

Fig. m 

The combination of overall and tangential compression 
leads to elongated structures that keep sweeping up cold 
gas. After ks 250kyr (Fig. [IJ middle panel) the dominant 
structures are already excavated by the combination of 
radiation and the pressure of the hot gas. From now on 
the evolution is mainly dominated by the hydrodynamic 
interactions between the hot and cold phase of the gas. 

After ks 500kyr (Fig. [TJ bottom panel) the morphol- 
ogy is remarkably reminiscent of the observed structures. 
The pillars in our simulations are indeed very complex 
structures with a cork-screw type, torqued morphology 
and show ro tational motion aro und their main axis, as it 
is observed (jGahm et al.ll2006h . Up to now it has been 
suggested that these complex morphologies arise due to 
magnetic fields, which are not included in our simula- 
tions. It is very likely that the pillars in M16 are a snap- 
shot of the formation scenario proposed here. At this 
stage the densest region (indicated by the center of the 
white box in Fig. [IJ bottom panel) undergoes gravita- 
tional collapse, the simulation is slowed down consider- 
ably and we terminate it. Future simulations with the 
inclusion of e.g. sink particles to avoid the detailed cal- 
culation of the further gravitational collapse leading to 
low mass stars will allow us to trace the subsequent evo- 
lution of the whole region. We call the most prominent 
feature in the white box in Fig. [T] (bottom panel) 'pil- 
lar I' and the second largest 'pillar IF, the collapsing 
compact core is at the tip of pillar II. Their respective 
masses are M pi n ar i = 12.3M©, M p m ar n = 8.1M© and 
M core = 0.7M Q . The compact core is defined as all mate- 
rial with a number density above n C rit = 10 7 cm~ 3 in a re- 
gion of i? C rit = 0.02pc around the peak density (see G08). 
Observations show that star formation is taking place 
close to the tips of the evolving structures (jSnider et al.1 
2007). The same is true for our simulations. In the pro- 
cess of sweeping up the dense material lags behind, gain- 
ing less momentum and thus leading to very high density 
enhancements near the radiation front. In contrast, the 
simulation without UV-radiation does not show any signs 

1 This timescale is calculated by taking the sound speed of the 
hot ionized gas c s h ot = 13.1pc/Myr and the average penetration 
length of the ionization of 0.4pc into account. 




l°gio k 



Fig. 2. — Evolution of the density weighted power spectra, plot- 
ted is the logarithm of the wavenumber versus the logarithm of the 
specific kinetic energy. Solid lines denote the run with ionization, 
dotted lines the control run without ionization. Plotted is the mean 
value of the spectra in the four different backward cubes, the error 
bars show the minimum and maximum values. Blue: total specific 
kinetic energy, green: solenoidal modes, red: compressional modes, 
dashed: Kolmogorov power-law. 

of gravitational collapse. 

Overall the scenario is very similar to the 'collect and 
collapse' model, as the denser regions would not collapse 
on their own on the timescale simulated and the sweeping 
up of material plays a vital role. We call it 'collect and 
collapse with turbulent seeds'. 

3.2. Turbulent Evolution 

For the discussion of the evolution of the density 
weighted spectra from the run with ionization we per- 
form a control run with the same initial conditions and 
accuracy parameters as the main simulation but with- 
out the inclusion of UV-radiation. In the comparison 
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run no sign of star formation can be detected. This is 
reasonable, since the cloud is not set up to become grav- 
itationally unstable as the total time of the simulation is 
much less than tg. The comparison is shown at t — Okyr, 
t = 250kyr and t = 500kyr (Fig. HJ. We analyze the tur- 
bulent spectrum in four cubic domains in the backward 
domain of the simulation, spanning 2pc in each direction 
(one is indicated in Fig[T] top panel, the other three are 
shifted in the negative y- and z-direction, respectively). 
Thus, it is guaranteed that there is always enough cold 
gas in the volume to achieve sensible results. To avoid a 
bias by either the ionized gas or the forming high density 
regions we take into account only gas with a number den- 
sity 10 2 cm~ 3 < n < 10 4 cm~ 3 . The particles are binned 
on a (128) 3 cubic g rid by using a kernel- w eighted binning 
routine (as e.g. in iKitsionas et al.ll2008f ). Based on this 
grid we calculate the density weighted spectra by sub- 
stituting v with (p/p) 1 / 2 (v — vrms) before the Fourier 
transformation, vrms is the average velocity in each of 
the three components. The specific turbulent kinetic en- 
ergy in the Fourier space is then given as 

£kin,turb = ' V '): ( 2 ) 

where v' and v' are the Fourier transform of the sub- 
stituted velocity and its complex conjugate, respectively. 
By mapping the ekin.turb cube to wave numbers k, the 
specific energy in the compressional, curl-free modes can 
be calculated as 

(v'-k)(v'-k) 

fc com — fc turb,kin . w , ■ \y) 

(k • k)(v' • V'J 

The specific energy in the solenoidal or incompress- 
ible, divergence-free modes is then given by e so i = 
ekin,turb — e C om- We construct the spectra by collect- 
ing the energy in the different wavenumber intervals (see 
IKitsionas et al.l [20081 for details). The total spectra as 
well as the solenoidal and compressional parts are shown 
in Fig. 1 

The initial spectrum at t = Okyr (Figj2l top panel) 
resembles quite well a power-law, even though the large- 
scale (low k) modes are lower, as the initial conditions 
are not produced by driven but by freely decaying turbu- 
lence. The slope is similar to a Kolmogorov-law. Approx- 
imately 25% of the total turbulent energy is contained in 
the compressional modes. 

At t = 250kyr (Fig[2j middle panel) there is already a 
distinct difference between the two spectra. The control 
run keeps the power-law shape and dissipates energy. In 
the ionization case the power is strongly increased on all 
scales. The increase is pronounced for k > 10, corre- 
sponding to scales < 0.2pc. An interesting feature is the 
rise in the compressional modes, where now 39% of the 
total turbulent energy is contained, whereas in the com- 
parison run this ratio stays at 25%. This clearly shows 
that the energy of the radiation is transferred into com- 
pression of the cold gas via the hot gas. The increase is 
in the turbulent energy itself and not correlated to the 
overall bulk motion in the x-direction, since the mean 
velocity is subtracted separately in each direction before 
the Fourier transform. 

After t — 500kyr (FigEJ, bottom panel) these differ- 
ences are even more pronounced. The kinetic energy in 
the cold gas is now a factor of four higher than in the 



run without ionization. Approximately 33% of the total 
turbulent energy is now contained in the compressional 
modes. This suggest that after an initial phase of high 
compression the system starts to relax. 

Including the mass in the respective region and den- 
sity range the total turbulent energy can be calculated. 
The initial turbulent energy is Bturb = 2.1 x 10 45 erg, 
the final turbulent energy (at t = 500kyr) is E[ on — 
4.3 x 10 45 erg and E n \ on = 1.1 x 10 45 erg in the ionized 
and unionized case, respectively. Thus, the input of tur- 
bulent energy per unit volume and unit time averaged 
over the simulation time when comparing the run with 
ionization to the case of freely decaying turbulence is 
etmb = 2.1 x 10 _25 ergs _1 cm -3 . By using the simplified 
assumption that the UV-radiation is absorbed isotropi- 
cally in the entire simulation volume the amount of en- 
ergy contained in the ionizing radiation for the chosen 
flux Fl y is eLy = 3.5 x 10~ 20 er gs~ 1 cm~ 3 . Compared to 
the es timates of lMatznerl (|2002| ) and lMac Low fc Klessenl 
(pool our radiative energy is several orders of mag- 
nitude higher, since we look at the direct surround- 
ing of an O-star instead of averaging over an entire 
galaxy. Nevertheless, the conversion efficiency of ioniza- 
tion into turbulent motion of the cold gas is in our case 
c = eturb/eLy ~ 2 x 10 -5 , which is an order of magnitude 
higher than their estimate of a « 2 x 10~ 6 for the Milky 
Way. Our highly resolved simulations show that ionizing 
radiation from an O-star or association provides a much 
more efficient mechanism to drive and sustain turbulence 
in the parental molecular cloud than was previously es- 
timated. However, this is still the energy input into the 
local environment in contra st to the average input rat e 
on galactic scales derived bv lMac Low fc Klessenl (|2004f ). 
On the larger scales it does not appear to contribute as 
significantly as e.g. supernova explosions. 

4. DISCUSSION 

We have shown in this letter that the observed pillar- 
like structures around O-stars as well as the gravitational 
collapse at the tip of the pillars can result from the im- 
pact of the ionizing radiation of massive O-stars on a 
turbulent molecular cloud. In addition, the turbulent 
energy in the cold gas is increased by a factor of four, es- 
pecially in the compressional modes. Both effects are due 
to the same mechanism: the ionization can heat the gas 
along channels of low density, thereby compressing gas 
at higher density into filaments. Close to the source of 
ionization this leads to the excavation of pillar-like struc- 
tures with triggered gravitational collapse at their tips. 
Further away from the source front, the structures have 
not yet fully developed, nevertheless the effect of com- 
pression is clearly visible in the turbulent energy spectra. 

Even though we find striking similarities between our 
simulations and observations, one has to bear in mind 
that this is a simplified approach which does not involve 
full radiative transfer. Ionized gas which gets shaded is 
assumed to cool immediately without affecting the ad- 
jacent structures. In addition, the shaded gas does not 
get ionized and heated by the recombination radiation 
of the ionized gas surrounding it. This might influence 
the precise shape of the structure behind the tip. More- 
over, the thin surface layers around each pillar where 
cold and hot gas are mixing cannot be resolved, although 
they might be crucial for the precise understanding of the 
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temperature and the chemical composition of these struc- 
tures. Nevertheless, our simulations indicate that these 
detailed effects arc of minor importance to explain the 
global picture, i.e. the overall structure and mass assem- 
bly of the pillars observed. Stellar winds might have an 
additional impact. Although O-stars have very powerful 
winds which can reach velocities of up to lOOOkm/s, our 
models suggest that ionizing radiation alone can repro- 
duce most observed features. 

The straightforward combination of hydrodynamics 
and ionizing radiation together with a standard turbulent 
model and typical parameters for molecular clouds leads 
to morphologies consistent with observed objects like pil- 
lars and collapsing cores. The similarities suggest that 



ionizing radiation plays a major role not only in shaping 
the parental cloud, but also in triggering secondary star 
formation. Furthermore, the overall turbulent kinetic en- 
ergy in the cold gas is increased significantly. 
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